Anomalous Interlayer Exciton Diffusion in WS2/WSe2 Moiré Heterostructure

Stacking van der Waals crystals allows for the on-demand creation of a periodic potential landscape to tailor the transport of quasiparticle excitations. We investigate the diffusion of photoexcited electron–hole pairs, or excitons, at the interface of WS2/WSe2 van der Waals heterostructure over a wide range of temperatures. We observe the appearance of distinct interlayer excitons for parallel and antiparallel stacking and track their diffusion through spatially and temporally resolved photoluminescence spectroscopy from 30 to 250 K. While the measured exciton diffusivity decreases with temperature, it surprisingly plateaus below 90 K. Our observations cannot be explained by classical models like hopping in the moiré potential. A combination of ab initio theory and molecular dynamics simulations suggests that low-energy phonons arising from the mismatched lattices of moiré heterostructures, also known as phasons, play a key role in describing and understanding this anomalous behavior of exciton diffusion. Our observations indicate that the moiré potential landscape is dynamic down to very low temperatures and that the phason modes can enable efficient transport of energy in the form of excitons.

Van der Waals (vdW) heterostructures are often considered key to unlock the full potential of two-dimensional (2D) materials. 1,2Weak interaction at the interface makes it possible to virtually stack any 2D system on top of another, creating artificial heterostructures with entirely novel properties.Twisted bilayer graphene (TBG), for example, displays flat bands near the Fermi level. 3The heavy Fermions populating these states not only give rise to unconventional superconductivity 4 but also insulating behavior driven by electron correlations. 5Strong electron correlation is also observed in transition metal dichalcogenide (TMD) heterostructures like WS 2 /WSe 2 , where flat bands arise 6,7 and a Mott insulator state and generalized Wigner crystallization are observed. 8,9These and many other emergent properties, 10−12 typically unexpected in the constituent monolayers, have been attributed to the moirépotential created by the superposition of two layers with different lattice constants and/or relative crystal orientation. 13,14In this paper, we show that the moirépotential is not stationary, and its dynamical nature is inferred by employing the temperature-dependent propagation of interlayer excitons (IX) in WS 2 /WSe 2 as the probe for the moirép otential dynamics.Our experiments suggest that the moireṕ otential is dynamic at finite temperature and the quasiparticle describing its motion, the phason, 15,16 strongly affects the exciton diffusion mechanism at low temperatures.
The WS 2 /WSe 2 system is made up of stacked direct band gap monolayers and can host both intralayer excitons and IX. 17,18The latter is formed when the electron and hole reside in different layers.The IX's lifetime is extremely long 19,20 compared to the few ps measured for monolayer TMDs 21 and can reach hundreds of nanoseconds in some systems. 22Its spectroscopic and dynamical features are dramatically affected by the presence of a moirépotential, as opposed to the exciton diffusion in a single monolayer where such moirémodulation is not present. 23,24An initial study via transient absorption spectroscopy of IX dynamics and diffusion in CVD-grown WS 2 /WSe 2 suggested that the moirépotential can trap excitons inhibiting their diffusion. 25Furthermore, a trapping/detrapping process, which is a temperature-activated mechanism where thermal energy allows the exciton to hop across the moirépotential wells, has been proposed for a MoSe 2 /WSe 2 heterostructure at various twist angles. 26Our observations of IX diffusion cannot be explained by such a model alone.For that reason, we consider the role of phonons, specifically of moire-induced phonons also known as phasons, 15 that are still experimentally unexplored.Recent calculations demonstrate how these quasiparticles can play a crucial role in understanding the charge ordering and correlated states of magicangle TBG. 27,28e perform a wide-ranging temperature-dependent study of the IX diffusion as observed through photoluminescence (PL).We explore the PL spectra of the interlayer excitons as a function of temperature and stacking angle (0°and 60°), together with the lifetime and diffusion length from which the diffusion coefficient is derived. 29This allows us to gain insight, first, into the temperature-dependent spectral properties of the IX for the two stacking configurations with different moirél andscape and then into the IX propagation dynamics.Our findings go beyond the classical Arrhenius model of thermally activated exciton transport involving a simple fixed potential landscape that have been reported in previous works where exciton diffusion is negligible at low temperatures (below 100 K). 25,26,30 Instead, we find that below 100 K excitons can explore the system, despite being fully trapped inside deep moirépotentials.An accurate model describing the moireṕ otential landscape is also provided, matching the energy barrier obtained from temperature-dependent diffusion measurements.
Finally, a mechanism that describes the electronic properties and vibrational modes of the two systems is proposed, showing how the moirépotential itself is affected by temperature and is likely the main mechanism for a nonzero diffusion coefficient at cryogenic temperatures, where the excitons are fully trapped.

RESULTS & DISCUSSION
The IX is an electron−hole pair where the two charged particles are in different layers, and the recombination is affected by the selection rules dictated by the twist angle.In fact, the different twist angles affect the electronic properties of the heterostructure, dictating the alignments of the spin-valley locked bands.To visualize the IX spectral properties, we analyze the PL emission as a function of temperature.Figure 1a,b reports the PL spectra of the two WS 2 /WSe 2 heterostructures with 0°and the 60°alignment in blue and orange, respectively (Figure S1).The PL spectra are reported at increasing temperatures, starting from 30 K. For both stackings, the emission is entirely dominated by the IX involving the top of the WSe 2 valence band (hole) and the bottom of the WS 2 conduction band (electron) at K and K'.Such emission is centered at 1.42 eV for the 0°stacking and 1.44 eV for the 60°stacking at 30 K, and progressively, red shifts with increasing temperature, as expected. 31,32It is important to note that the 0°and 60°stacking configurations exhibit different dielectric screening due to the distinct electronic environments that they create.Consequently, there is a variation in their exciton binding energies, 33,34 highlighting the need for careful comparison between these structures.Moreover, the exfoliation and stacking method employed in device fabrication also introduces variability in coupling strength between layers , potentially leading to observable differences in emission energy.Interestingly, for the 0°s tructure, two thermally activated IXs emerge: a first one at around 75 K and a second one above 150 K.The insets of Figure 1a,b show representative spectra featuring all peaks at three selected temperatures.We can describe the peak evolution by tracking their position as a function of temperature, as displayed in Figure 1c.We evaluate each peak position by fitting the spectra with one, two, or three Voigt functions, depending on the number of observed peaks (Figure S2).Using Varshni's equation, which empirically describes the temperature dependence of energy gaps in semiconductors, 32 we extract an empirical electron−phonon coupling coefficient associated with the average phonon energy ⟨ℏω⟩.The fitted parameter values are reported in the inset of panel (c).The two lowest excitons from the 0°structure (marked with blue circles and triangles in Figure 1a) have similar fitting parameters, suggesting very similar electron− phonon coupling.Moreover, they display an energy separation of ΔE ∼ 30 meV, compatible with the spin orbit splitting of the conduction band in WS 2. 35 The origin of the third, highest energy peak marked by a blue cross, which displays a much steeper slope, is unclear as of now; however, it could be associated with one of the moiréexcitons previously observed in the same kind of heterostructure for the intralayer exciton. 11nother plausible hypothesis considers the description of the moirépotential as a quantum well; within this hypothesis, the third peak arises from the optical transition involving the first excited state of the well, similarly to what reported by Chatterjee and co-workers. 36The PL spectra collected for the 60°stacking sample show different behavior.There is only one peak (marked with red diamond), and it shows an identical trend as the two lowest energy emission peaks from the 0°sample.The presence of three peaks in the 0°alignment and the single peak emission for the 60°can be ascribed to the selection rules coming from the electronic structure of the single layers. 37,38The lowest energetic transition for the 0°s ample should be related to the spin-dark exciton.However, the presence of the moirépotential has been demonstrated to partially lift those selection rules in systems with large lattice mismatch, 37−40 and given it is the lowest energy state, its high population at low temperatures makes the PL emission comparably strong, similarly to what observed for MoSe 2 / WSe 2 structure. 41With increasing temperature, however, the tail of the Fermi−Dirac distribution starts populating the higher spin-split level of the WS 2 conduction band, allowing the spin-aligned IXs to become observable. 18On the other hand, the 60°sample has a bright transition as the lowest available one, allowing all the hot electrons to relax and recombine through that channel, while the higher spin-dark state, being much less efficient and much less populated even for elevated temperatures, is not easily observed (Figure 1d).It is crucial to note that the extended lifetime of excitons, as discussed later in the paper, further facilitates the relaxation from the spin-bright to the spin-dark state in the 0°sample.This longevity allows excitons more time to undergo this transition, augmenting the population of the spin-dark state and influencing the observed PL emission patterns.
We now focus our attention on IX dynamics, studying its diffusion and recombination lifetime.We perform an exciton diffusion study by collecting the light associated with the IX transitions only and mapping it in real space (Figures S3 and  S4).We use a 790 nm long pass filter in detection to cut out any potential light emitted from the intralayer transitions, particularly at elevated temperatures.
The diffusion length is defined as 42 where σ L is the Gaussian covariance that fits the centrosymmetric excitation laser profile along an arbitrary axis (Figure 2a) and σ is the Gaussian covariance of the real space PL intensity, evaluated in the same manner, assuming that, if excitons can diffuse, the emission of the initially excited 2D Gaussian profile broadens isotropically over time, following Fick's second law.The diffusion length L is evaluated as the square root of the difference between the squared covariances, according to eq 1 (Figure 2b).Simultaneously, we collected the time-resolved PL emission.Figure 2c,d shows data for the 0°and 60°samples, respectively.The decay traces exhibit a single exponential trend, which we fit to extract PL lifetime τ (Figure 2e).
The diffusion length for both stacking angles is largely unaffected by temperature variations.It is mostly confined within a mean value of ∼600 nm, which is consistent with reported diffusion lengths in monolayer systems. 43,44On the other hand, the PL lifetime shows a clear decreasing trend with increasing temperature, dropping after 150 K for 0°and 100 K for 60°.Using these values of L and τ, we extract the diffusion coefficient, defined as 43,45 Figure 2f shows the calculated values of diffusion coefficient D (diffusivity), which undergoes a thermally activated process resulting in a pronounced increase of the diffusivity with increasing temperature.Such an exponential behavior is expected when considering a simple moire-trap picture. 25,26,30nce formed, the IXs are affected by the presence of the moireṕ otential landscape and are trapped in its local minima.However, if enough thermal energy is provided and the activation barrier is overcome, the IX can diffuse faster, reaching a value comparable with the one of intralayer excitons of an isolated monolayer where the moirépotential is not present (Figure S5).To evaluate the effective potential for diffusion, we perform extensive density functional theory (DFT) calculations for both 0°and 60°twisted heterobilayers, considering the atomic reconstructions at T = 0 (see the Supporting Information).The moirépotential is evaluated as E g (r) represents the local bandgap and E g min represents the averaged bandgap. 37Since Δ M (r) varies smoothly over the moiréperiod, the effective moirépotential for diffusion can be extracted from the local bandgap at different high-symmetry stackings.The averaged effective moirépotential for 0°(124 meV) is unmistakably larger than that of 60°WS 2 /WSe 2 (85 meV) (Figure 2g,h).
The most striking observation, however, is that a simple detrapping model is not consistent with what observed at low temperatures (T < 100 K). Figure 2i displays the diffusivity data reported in the dashed box in Figure 2f.We emphasize that an Arrhenius-like description would make D tend to 0 at low temperatures (dashed line).However, the data clearly show a finite offset D 0 , indicating that the IXs are diffusing, even at low temperatures.The diffusivity as a function of temperature now takes the form: The moirépotential height E 0 can therefore be extracted from the data resulting in a value of 126 and 67 meV for 0°and 60°orientation, respectively, very similar to the theoretically calculated potential reported above and shown in Figure 2g,h.Values of D 0 = 0.16 and 0.24 cm 2 /s are also obtained from the fitting for 0°and 60°, respectively.To understand the nature of D 0 , we aim to juxtapose our findings against existing models detailed in the literature.−48 These studies have contributed to understanding exciton transport mechanisms within these materials, including effects such as weak localization.Our work extends this exploration to heterobilayer TMDs, a domain where the moirépotential, not present in monolayer systems, notably influences IX behavior.In our heterobilayer system, exciton transport at high temperatures is characterized by thermally activated hopping across moirépotential minima, demonstrating an Arrheniustype temperature dependence.This established mechanism, however, does not account for the observed finite diffusivity at lower temperatures.Unlike monolayers where phonon scattering could facilitate exciton transport, the heterobilayer's thermal phonon energy scales fall significantly below the moireb arriers evaluated above, suggesting an alternative transport pathway.Other effects such as impurities, defects, and disorder play a less critical role in our heterobilayer observations.This conclusion is supported by other analyses showing that these elements would typically lead to exciton trapping 49 and reduced diffusivity at low temperatures contrary to our findings of maintained finite diffusivity, suggesting minimal trapping within our system.Furthermore, our examination of quantum effects, supported by low exciton density estimations, indicates the absence of significant exciton−exciton interactions.In fact, we estimate the exciton density (as reported in the Experimental Methods section) to be significantly low (5 × 10 10 cm −2 ), much smaller than one exciton per moiréunit cell.This assessment is corroborated by the photoluminescence spectrum of the WSe 2 monolayer, which lacks any signs of exciton−exciton peaks 50 under experimental conditions similar to those of our IX studies (see Figure S6).Other effects such as exciton−electron interaction that have been demonstrated to lead to an increase in the diffusion coefficient with increasing carrier density 48 can also be excluded.In fact, the PL spectrum of monolayer WSe 2 displayed in Figures S5 and  S6 does not display the typical feature of the charged exciton peak arising from the presence of high free carrier density. 51e propose moiréphonons (phasons) 15 to be responsible for the deviation from the Arrhenius description.The phason modes represent an effective translation of the moirésites due to the nonuniform out-of-phase in-plane translation of two constituent layers in the heterobilayer.A mechanism described as "moirésurfing" has been proposed to affect electron transport in twisted TMD homobilayer. 52It is unique to moireś ystems, where thermally populated phason modes enable the collective motion of the moirélandscape, thereby, in our case, transporting excitons.This process diverges fundamentally from the scattering phenomena previously considered for monolayers, offering another perspective on exciton mobility in complex TMD systems.
Here, we calculate that the phason modes give rise to dynamic moirépotentials at finite temperature.We simulate a WS 2 /WSe 2 heterostructure in both twist angle configurations in a canonical ensemble using classical molecular dynamics (MD) simulations.We find that the moirépattern starts to move in-plane due to the time-dependent relative local displacement of the two layers, as a result of the thermally activated phasons.Figure 3a−h shows different temporal snapshots of the valence band maximum (VBM) and conduction band minimum (CBM) single particle density of states for 0°and 60°stacking at a lattice temperature of 175 K.The MD simulation reveals, as evidenced by the colored marks, that the different stacking sites, and thus the moireṕ otential, move with time, affecting the exciton hopping (Figure 3i).It is also shown that the single particle density of states calculated from DFT strictly follows this movement as would be expected within the Born−Oppenheimer approximation (see Supporting Information Figure S9.An animated simulation is found in the online Supporting Information representing the crystal dynamics over 2 ns at 175 K).
As a preliminary indication of our model's accuracy, we compare the speed of phasons derived from phason dispersion calculations (Figure 4a,b) and MD simulations with the experimentally observed value.By calculating the phason group velocity directly from the phonon dispersion curve (Figure S10), we estimate the speed to be up to 200 m/s.Note that the phason has a parabolic dispersion close to the Γ point of the BZ corresponding to a vanishing group velocity reflecting the small, but finite energy to slide a commensurate moirébilayer.It is also possible to extract the speed of the moiréphason as a function of temperature from the MD simulation by evaluating the stacking site displacement at different times (Figure 4c).The speed increases sharply with temperature and then quickly saturates, in agreement with the expectation that the sliding of two dissimilar layers costs little energy, which is generally referred to as super lubricant behavior. 53,54The value obtained from the MD simulations is considerably lower, with respect to the fitting of the phason acoustic branch, approximately 20 m/ s.The reason behind this mismatch is rooted in the parabolic dispersion around Γ, deviating from a typical linear phonon acoustic branch, whose slope is the speed of sound.It is notable that the speed estimated from the experimental measurements�obtained by dividing the diffusion length by the radiative lifetime�is found to be around 50 m/s, which lies within the bounds set by the two theoretical methodologies described above.The agreement of our experimental results with the boundaries defined by our simulations underscores the necessity of careful handling and interpretation of experimental results when dealing with moirésystems.

CONCLUSIONS
This paper studies the IX diffusion in WS 2 /WSe 2 heterostructures at two different stacking angles as a function of temperature.The distinct PL signatures we have observed are tied to the intrinsic selection rules governing the monolayers that constitute these heterostructures.Intriguingly, we have determined that the moirépotential appears to exhibit different depths for the two stacking angles.Specifically, a shallower potential of 67 meV is associated with the 60°stacking, while a deeper 126 meV potential characterizes the 0°stacking.Unexpectedly, the most striking observation is the existence of an anomalous nonzero diffusion at cryogenic temperatures in samples of both orientations.This cannot be explained within the conventional framework of an activated process that detraps the excitons from the moirépotential well.To make sense of this unconventional behavior, we suggest that IX may be effectively transported across the crystal lattice by a moving moirépotential landscape.Our MD simulations indicate that this process may initiate from temperatures as low as 5 K. Lower temperature studies of IX dynamics may further elucidate the complex relationship between excitons and the phonons of moirésystems.The consistency we have observed between theory and experiment, considering the matching moirépotential barrier and moiréspeed, serves as a compelling indication that moiréphonons, phasons, are responsible for the anomalous diffusion observed in this work.These findings underscore the critical significance of incorporating phason dynamics into forthcoming investigations of moirésystems.

EXPERIMENTAL METHODS
The heterostructure is prepared with a standard pick up technique. 55The TMD heterostructure is incapsulated in h-BN to preserve the pristine properties of the system.The relative alignment has been checked using second harmonic generation (Figure S1).The experiment has been carried out using a Montana Cryostation s200 closed-cycle helium cryostat with full temperature control down to 4.2 K and 100x optical magnification (NA = 0.95).The samples are excited using a pulsed (38.9 MHz) supercontinuum laser tuned to 510 nm with 10 nm line width and 30 ps pulse length, focused down to a spot diameter of about 320 nm fwhm, yielding an energy density per pulse of 0.38 μJ/cm 2 for our diffusion experiment.Assuming 5% absorption mainly from the WS 2 B-exciton and unity conversion into interlayer excitons, this amounts to an exciton density of roughly 5 × 10 10 cm −2 , within a regime where we can exclude a nonlinear response.At this density, quantum effects due to dipole−dipole interaction are negligibly small compared to the thermal energy in the range of temperatures explored. 56PL spectra are collected using a spectrometer with a dispersion grating together with a Peltiercooled charge-coupled device (CCD) camera, while lifetimes are measured using an avalanche photo diode (APD) with a nominal temporal resolution of 30 ps.Real space imaging of the PL emission to extract exciton diffusion lengths is realized using a fast and highly efficient CMOS detector in the time integrated detection mode.

COMPUTATIONAL METHODS
Twisted WS 2 /WSe 2 heterobilayers have been generated using the TWISTER package. 57All of the structural relaxations, molecular dynamics, and phonon calculations are performed using classical interatomic potentials fitted to DFT calculations.The intralayer interactions within WSe 2 and WS 2 are described using Stillinger-Weber potential, 58 and the interlayer interactions are captured using the Kolmogorov-Crespi potential. 58All the simulations with classical interatomic potentials are performed using the LAMMPS package, 59 and the phonon calculations are performed using a modified version of the PHONOPY package. 60We relax the atoms within a fixed simulation box with the force tolerance of 10 −6 eV/Å for any atom along any direction.All the MD simulations are performed in the canonical ensembles for several temperatures, and the moirémovements are extracted in the microcanonical ensemble.The speed of moirésites is extracted from MD simulations of a 3 × 3 × 1 supercell of the moiréunit cell.We find that the moirésites move even for a simulation of a 20 × 20 × 1 supercell.Note that the supercell size dictates the lowest accessible momentum close to Γ point, i.e., q min ∝ 1 L m , where L m is the moirésimulation cell lattice constant.The moirésite speed differences between twist angles 0°and 60°a re reproduced for 1 × 1 × 1 and 6 × 6 x 1 supercells as well.All the electronic structure calculations are performed using the SIESTA package. 61All the calculations include spin−orbit coupling. 62We have used the local density approximation with the Perdew−Zunger parametrization as the exchange-correlation functional. 63

ASSOCIATED CONTENT Data Availability Statement
The twisted heterobilayer structure construction, atomic relaxations, molecular dynamics simulations, phonon spectra calculations, and electronic band structure calculations presented in the paper were carried out using publicly available codes.Our findings can be fully reproduced using these codes.Some of the postprocessing tools, such as tracking the motion of moirésites, developed for the paper will be made publicly available at https://gitlab.com/_imaity_/Moiredynamics.search Square.https://www.researchsquare.com/article/rs-2627775/v1(accessed May 9, 2024).

Figure 1 .
Figure 1.PL spectra of 0°(a) and 60°(b) stacking angles at different temperatures.Insets display representative spectra collected at different temperatures.Blue and red symbols identify the different peaks present in the spectra.In panel (a), a WSe 2 intralayer exciton peak is highlighted at higher temperatures.(c) Peak positions tracked as a function of temperature, fitted with the Varshni equation (black line).Inset displays the value of the average phonon energy for each fit.(d) Scheme of the spin-polarized band structure for the 0°and 60°(left and right panels, respectively) structures with the corresponding vertical optical transitions (dashed lines).

Figure 2 .
Figure 2. (a) Exciton diffusion length evaluated from PL area with respect to the laser excitation width at 30K.(b) Diffusion length values at different temperatures for the two stacking angles.PL lifetime traces at different temperatures for 0°(c) and 60°(d) stacking angle.(e) PL recombination lifetime as a function of temperature for the two stacking angles.(f) Diffusivity values fitted with (solid line) and without (dashed line) D 0 offset.(g) and (h) moirépotential evaluated as described in the main text for 0°and 60°, respectively.Dashed black diamond indicates the moiréunit cell of periodicity λ ∼ 7 nm.(i) Zoom-in of the diffusivity map at a lower temperature range to highlight the D 0 .

Figure 3 .
Figure 3. Local holes and electron densities evaluated at valence band maxima (VBM) and conduction band minima (CBM), respectively, for 0°and 60°deg structure.VBM and CBM densities are evaluated at different times, t 1 (a−d) and t 2 (e−h) for finite temperature T = 175 K with t 2 − t 1 = 0.3 ns.Colored dots and diamonds track the different moirésites of the heterostructure.(i) Scheme of the exciton (yellow sphere) transport (blue arrows) assisted by moiréphasons (red arrows).

Figure 4 .
Figure 4. (a,b) Phason dispersion overlaid to the phonon dispersion for 0°and 60°, respectively, computed at T = 0 K for the longwavelength phonons near the Γ point.Phason speed is reduced by a factor of 2 and 1.6 when compared to the LA mode for 0°and 60°twist angles, respectively.(c) Speed of moirésites captured using molecular dynamics simulations as a function of temperature for 0°and 60°.